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We investigate the long time behavior of a wavepacket initially localized at a single site no in 
translationally invariant harm onic an d anharmonic chains with random interactions. In the har- 
monic case, the energy profile {e n (t)) averaged on time and disorder decays for large |n — no| as a 
power law {e n (t)) ~ C\n — no\~ v where n = 5/2 and 3/2 for initial displacement and momentum 
excitations, respectively. The prefactor C depends on the probability distribution of the harmonic 
coupling constants and diverges in the limit of weak disorder. As a consequence, the moments 
(m v {t)) of the energy distribution averaged with respect to disorder diverge in time as t^'"' for 
v > 2, where /3 = v + 1 — r\ for v > r\ — 1. Molecular dynamics simulations yield good agreement 
with these theoretical predictions. Therefore, in this system, the second moment of the wavepacket 
diverges as a function of time despite the wavepacket is not spreading. Thus, this only criteria often 
considered earlier as proving the spreading of a wave packet, cannot be considered as sufficient in 
any model. The anharmonic case is investigated numerically. It is found for intermediate disorder, 
that the tail of the energy profile becomes very close to those of the harmonic case. For weak and 
strong disorder, our results suggest that the crossover to the harmonic behavior occurs at much 
larger |n — no| and larger time. 

PACS numbers: 05.45.-a 05.60.-k 42.25. Dd 



I. INTRODUCTION 



There has been large activity for many years in the study of the temporal evolution of an initially localized energy 
excitation in various nonlinear systems, e.g. the discrete, nonlinear Schrodinger equation (DNLS) [lj-|4j, Fermi-Pasta- 
Ulam (FPU) and Klein-Gordon (KG) model |, || with both uniform and random couplings. In the latter case, 
the main interest is in the interplay of anharmonicity (nonlinearity) and disorder which is not yet fully understood. 
For harmonic one-dimensional disordered systems, all eigenmodes (called Anderson modes) of the infinite system are 
known to be localized and form a complete basis. Then a wave packet at time t = will remain localized at any time 
as a linear superposition of Anderson modes of the infinite chain. Whether or not this behavior changes qualitatively 
by introduction of anharmonicity is highly debated and controversial (see Refs. 0,0; an d references therein). 

Since an analytical treatment of the time evolution of anharmonic systems with disorder is extremely difficult, most 
investigations have been done by molecular dynamics simulations. In the numerical studies, one typically follows the 
wavepacket dynamics by monitoring quantities like the participation ratio P(t) (a measure the localization at time t), 
and the time-dependent moments m v {t) of the local energy e n (t) (see the definitions below). All this measurements 
are hampered by statistical errors as well as finite size and finite time effects. Even very long calculation times of, 
say, 10 s microscopic time units (of order picoseconds) may not be entirely conclusive. Indeed, one can never be sure 
whether the spreading of a wavepacket is complete or only partial in the infinite-time limit. This issues are intimately 
related to the spontaneous self-trapping of energy (for example in the form of discrete breathers) which is generic in 
most nonlinear systems. 

Independently on complete or incomplete spreading, one might expect that the evolution of the wavepacket tails 
should yield relevant information on the spreading process itself. In such regions, the typical displacement becomes 
small enough such that linear approximation of the forces becomes valid. This motivates the investigation of the 
harmonic chain, as a first necessary step for an insight of the nonlinear case. Despite the apparent simplicity of such 
a case there are still issues that have not been fully discussed in the literature. Let us briefly review some of the main 
results known for this case. Without disorder all eigenstates are extended and it is well-known (see e.g. Ref. [lfj, and 
references therein) that 



m u {t) ~t"M 



(1) 
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with (3(2) — 2, i.e. the energy spreading is ballistic (note that v is not necessarily an integer). Introducing disorder 
and/or anharmonicities, this energy transport is changed and may be superdiffusive ((3(2) > 1), diffusive ((3(2) = 1) or 
subdiffusive ((3(2) < 1), or it could become logarithmic or disappear ((3(2) = 0). If the initial excitation is at site zero 
with amplitude i*o(0), then the disorder averaged propagator (u n (t)) is one of the basic quantities. Although (uo(t)) 
for t — > oo is known analytically for different classes of disorder [ll|, much less is known for h ^ 0. Approximating 
the Anderson modes by plane waves with exponentially decaying amplitudes, it has been shown in Ref. |6| that 



, , ^^ 1 r (M — C.t) 2 . 
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for |n| — > oo and t — > oo. Here, c is the sound velocity and £o a measure of the localization length. Eq.Q shows for 
t — > oo there are two humps which propagate ballistically at the sound velocity c, but with an amplitude which decays 
as 1/ \fi. Within its co- moving frame , these humps spread as for normal diffusion. Another approach for calculating 
(u n (t)) is to use a scaling hypothesis (llj 



(5 n (u)) = (fi (w))f(n/(H) , uj^O (3) 

for the Laplace transform of (u n (t))) for ui —> 0. A similar Ansatz can be made for (u n (t)) Here, denotes 

a localization length. 

In this paper we investigate the energy profile (e n (t)) averaged on time and disorder of a wavepacket originating 
an initially localized excitation. We demonstrate that it asymptotically decays as a power law in space. Thus, the 
wavepacket remains localized only weakly while its moments appear to diverge in time. This result, which, to the 
best of our knowledge, has not been reported previously, must be taken into account expecially when attacking more 
difficult nonlinear case. Indead, some numerical results for the anharmonic chain (a FPU model) will be critically 
analyzed on the basis of the results on the harmonic one. 

The outline of our paper is as follows. In Section II we will introduce the harmonic model, rephrase some of its 
well-known properties, define the local energy e n (t) and give some information on our numerical approach. A virial 
theorem for the time averaged local kinetic and potential energy will be proven in Section III. It will be applied in this 
section for the analytical calculation of the time and disorder average of e n (t). The corresponding analytical result 
will be compared with the numerical one. Furthermore we will investigate the moments m u (t) of the local energy 
(s-n(t))- The influence of anharmonicity on e n (t) will be numerically studied in Section IV, and the final Section V 
contains a summary and some conclusions. 



II. THE DISORDERED HARMONIC CHAIN 



A. Property of the Anderson modes 

As motivated above we investigate the classical dynamics of a disordered harmonic chain with lattice constant a 
which is invariant under translations. Its classical Hamiltonian reads: 

(4) 

Here, u n is the displacement of the particle at site n, p n the corresponding conjugate momentum, m the particle's 
mass, and K n the random coupling constants between nearest neighbors. The K n are independent random variables, 
identically distributed with some probability distribution p(K). Stability requires all K n to be positive. In our 
numerical approach, the system is finite with N particles and with free ends, i.e. K± N / 2 = 0. Otherwise, we shall 
perform analytical calculations in the thermodynamic limit N — > oo where the choice of the boundary conditions does 
not matter. The equations of motion are 

mu n = K n (u n+ i - u n ) - K n -x(u n - u„_i). (5) 

The general solution of Eqs. ([5]) with initial conditions u n (0) = u n , u n (0) = u n is given by 



Pn 



7T~ + ^K n (u n+1 - u n ) 2 
2m 2 



u„(t) = U + U t + v n (t) 



(0) 
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where 

ri v=£Q i/^O v 

and Uq — U n/N, Uq = Y2 u n /N are the position and velocity, respectively, of the center of mass of the whole chain. 

n n 

The eigenmodes Qn with eigenfrequency ui v can be chosen as real with indices v in a countable set. They satisfy 

K n {QM Q%) + if„-i(QW - QWj = rn^QM ( 8 ) 

and they can be normalized, except the uniform eigenmode = 1 with loq = which is extended and cannot be 
normalized for the infinite system. For any size N of a finite system, the translation invariance of the model implies 
that Ql 0) = -7= is an eigenmode with eigenfrequency ujq = 0. In the limit of an infinite system, all eigenmodes 
are localized, except this single zero-frequency mode which is extended. However, nothing changes in the problem 
when choosing the center of mass of the whole system immobile at Uq = with Uq = 0. Though the eigenspectrum 
is discrete for the infinite system, it is dense. The corresponding density of states 



1 N ^ 

g{u) = lim — V 6(u - u v ) (9) 

N— foo iv ' — * 



is a smooth function which is known [15| to be self-averaging, i.e. it is independent on the disorder realization with 
probability one. Moreover, in the small frequency limit, u> — > 0, we have 0,113 



g(ui)~y-A / . (10) 

The localized eigenmodes, decay exponentially with a localization length [l|| [l7| 

which diverges at the lower "band" edge at ojq = 0. 

Then, if the chain is finite with length L, there is a frequency such that the localization length equals the system 
size, i.e £(wl) = L = aN. Consequently, only the eigenmodes with frequency ui v > lul can be considered as well 
localized inside the finite system, while the remaining modes where lu u < ojl extend over the whole finite system. 
Their number which is of order of \/N goes to infinity in the limit of an infinite system despite their relative weight 
for N — > 00 goes to zero as 1/y/W. As a result, they still play a role for transport quantities, like the energy diffusion 
constant 0, [T3, [H| or the thermal conductivity [l!| . Actually, those relatively extended modes behave like acoustic 
modes whose effective sound velocity is 



c=J<^a (12) 
V m 

Although these results were originally proven for a chain with mass disorder they also hold for our model. Indeed, 
letting y n = (it n +i — u n )/K n , Eq. ([5]) is mapped onto the eigenequation with mass disorder. This property has already 
been used above since the mass average (m) has been replaced by (K^ 1 ). 



We define the local energy: 



with kinetic and potential parts 



B. Local energy and local virial theorem 



3„(t) = e ( kin '(i)+ei pot) W 



e^ n) (t) = -(u n (t)) 2 (13) 
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and 

e { n pot) (t) = ^K n . x [u n {t) - u n ^{t)]u n {t) - 

-K n [u n+X {t) - u n (t)]u n (t) (14) 
respectively. Then, ei pot ' ) equals the total potential energy in Eq. @. We will investigate the energy profile for a 

n 

displacement excitation: 



u n {0) = AS n , no , p„(0) = (15) 



and a momentum excitation: 



u n (0) = , p n (0) = B6 nino . (16) 

For calculating numerically {u n (t)} we considered the example of a uniform and uncorrelatcd distribution of random 
couplings K n with probability distribution 

, x Crnf-FT Hk<K<Rk, 
p(K) = fc («-D ~ . ~ (17) 

I otherwise. 

where of course R>1. 

To explore the role of different disorder strengths we fixed k = 1 and took different R. The choice of units is such 
that tyl — 1 and a = 1. Note that with this particular choices the effective sound velocity, Eq. (JT2J) , is c = J jjjjj;- 

Microcanonical simulations were performed for typically A" = 8192 particles with fourth order symplectic algorithm 
[20j . with typical time step 5 x 10~ 3 or less. Although the choice of the initial conditions, Eqs. (fT"5"|) and (p~6|) . implies 
J7 = A/N, Uq — and Uq = 0, Uq = B/N, respectively, these nonzero quantities arc rather small, since A and B are 
of order one and N S> 1 . 

In our numerical experiments, we avoid that the wavepacket reaches the chain boundaries which may generate 
spurious finite-size effects (reflexions etc.). Thus, one should restrict the maximum simulation time t s - lm to be smaller 
than t m ax ~ N/ c where c is the sound velocity. 

We also fixed no = —N/2 + 1 for extending the spatial range of our system, so that one simulates the wavepacket 
propagation in a semi- infinite medium Q. Some runs with uq = where also performed, yielding similar results. 
Figure [1] shows the numerical profile e n (t = 2000) for a momentum excitation with B = 2.0. The result for a single 
realization of the disorder exhibits on the log-log-rcprcsentation strong fluctuations around an average, decaying 
linearly. Averaging over a large enough number [O(10 3 )] disorder realizations strongly reduces these fluctuations and 
supports the linear dependence on \n — no| on the log-log scale. In the next Section we demonstrate that this is indeed 
the case and compute analytically the exponents associated with such power-lay decay. 

The calculation of the time averaged energy profile will be simplified by means of a local virial theorem, that will 
be proved below. The well-known virial theorem [2l[ relates the time average of the total kinetic energy to the time 
average of the virial. The virial [2l| involves the gradient of the total potential energy. If the potential is harmonic 
this theorem implies equality between the time averaged total kinetic and total potential energy. In this subsection 
we will prove that this relationship also holds for the time averaged local kinetic and local potential energy, defined 
by Eqs. (flU]) and (|T4")) . respectively. 

The time average of a function f(t) is defined by 

T 

7W= r «m ^ / dtf(t) . (18) 
o 

Substitution of the general solution u n (t) of Eqs. ([6|), into Eq. ([13]) and taking into account 



COS LO^t COS UJ^rt = — b~ vv t 



sin w v t sin u> v it = — 5„„< 



sinu^t costJv't = . (19) 
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FIG. 1: (Color online) Energy profile e n (i) at t = 2000 for a momentum excitation with B = 2.0, N — 8192 particles, R — 4 
and for a single realization of disorder (green line) and averaged over 2 x 10 3 realizations (black line). 



yields 



(20) 



Note that the sum over v remains discrete and cannot be replaced by an integral in the limit of an infinite system. 
With our definition Eq. (|14[) of e„ , and using Eq. ([5]), we obtain: 



i pot) W = ~u n (t) Un (t) 



(21) 



Substituting u n (t) from Eqs. ([6]), ([7]), and since u n (t) and 'i'n(i) bas to remain bounded at all times for any initially 
localized wavepacket (with finite energy), yields 



,(p°t) 



(<) = el kin) (t) 



(22) 



for all n and arbitrary initial conditions with finite energy in case that the center of mass has been chosen immobile 
(U = 0). 



III. ENERGY PROFILE: HARMONIC CASE 



A. Energy profile 

Without restricting generality we choose m = 1 and a = 1. Let us discuss first the case of a displacement excitation 
for a given disorder realization. In this case, we obtain from Eqs ([6]) and {?]) for A = 1 and Uq = 0, Uq = 



un{t) = £)qMqM cos w„t 

u n {t) = -^u v Q^Q^wa.w v t (23) 

and therefore 

e£ kin) W = ^ E UMQPQttQPQW^tMU't ■ ( 24 ) 
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Let us discuss first the qualitative i-dependence of e„ (t) . We will explain how the spectral properties govern its 
time dependence. Particulary we show that this quantity which is not averaged over time and/or disorder does not 
decay for n — > oo and/or t — > oo. Since the eigenspectrum of the infinite random system is discrete with a countable 

basis of localized eigenstates {Q^}, u n (t) has been expanded in this basis (see Eq. ([7]) and Eq. (p?5j) '). This expansion 
is actually an absolutely convergent series of cosine functions of time because 

iEqWi < (E^ } 2 ) 1/2 (E^ 2 ) 1/2 = 1 

Consequently, u„(t) is an almost periodic function in the sense of H. Bohr [25| . An equivalent definition for such 
functions is that for any arbitrarily small e > 0, there is a monotone sequence of t p {p € Z) (called pseudopcriods) 
which is relatively dense (that is there exists L such that r p+1 — r p < L for any p) and such that for all p, f(t) is 
periodic with period t p at the accuracy e that is \f(t + t p ) — f(t)\ < e for any p and for all t. As a consequence of this 
recurrence property, an almost periodic function cannot go to zero for t — > ±oo. The set of almost periodic functions 
is an algebra, that is linear combinations and products of almost periodic functions are almost periodic functions, as 
well. 

In our case, the set of eigenfrequencies uj v is bounded (since the support of the distribution function p{K) is 
compact) and thus it is straightforward to show that the time derivative u n (t) is also an almost periodic function of 

time, and the local kinetic energy el in) (t) defined by Eq. (13]), as well. el kin) (i) fr om Eq. (|24[) can be decomposed 
into a time independent term and remaining time dependent terms : 



+- A Ew^(QMqM)(Q^Q^)[cob(w w -w v O* + ^K+w,,0<] • (25) 

Note that for a finite chain without disorder, i.e. K n = K, the first and second term on the r.h.s of Eq.([25|) are of 
order 1/N since the eigenmodes are plane waves where Q { n ] oc \/VN. Then {Q^Qnof oc 1/A 2 , and there are only 
N such terms. Consequently, they will not contribute to e^ m \t), in the limit N — > oo when the eigenspectrum of 

the chain becomes absolutely continuous. In that case e^ fcm ^ (t) can be represented by an integral which is a Fourier 
transform of a smooth function and is obviously not an almost periodic function. It decays to zero at infinite time 
as expected from ballistic diffusion. This is not true in case of disorder, because each term in the series keeps a non 
vanishing contribution for the infinite system and el kln ^ (t) does not decay to zero at infinite time because it is almost 
periodic. 

However, in contrast to the ordered chain, Qn^QriJ 1S not a smooth function of w„, in case of disorder. The 
reason is that when the eigenspectrum is discrete, arbitrarily small variations of lu u may change the location of the 
corresponding localized eigenstate by arbitrarily large distances. Thus, these eigenstates {Qn^} arc not continuous 

functions of uj v but depend on the disorder realization as well as ei km ^ (t) and e n (t) (since they are obtained as discrete 
series explicitly involving these eigenstates). The consequence is that those quantities are not self-averaging, as clearly 

demonstrated by Fig. [1] for e„ (t). 



B. Disorder averaged profile 

Since e[i %n ^ (t) is an almost periodic function of time, it is a stationary solution. Its time average drops all cosine 
terms in Eq. (|25[) and keeps only the constant term, i.e. we get 



An attempt to justify the use of the time averaged quantity will be given below, ej, (i) and e n (t) still depend on the 
disorder realization. Therefore it is reasonable to calculate the corresponding disorder averaged quantities, as well. 
Despite they cannot be observed for any single disorder realization, they give information on the general behavior of 
the profiles. Then we arrive at 
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(ef n) (t)) = ^E^(^ ) ^o ) ) 2 ) (26) 

for the infinite system. 

Note that, in the infinite system, the set of eigenvalues and eigenvectors are discontinuous functions of the disorder 
realization. Yet, according to Wegner 22|, the disorder average (F{{Qn })) of an arbitrary function F({Qn)}) of 



the eigenvectors can be well defined as a smooth function of uj as a limit for finite systems with size N — > oo 



N 



\ w<w,y<w+5u; / n=l 



The sum in the integral is restricted to eigenvalues ui v which belong to an interval [w,ui + Slo] of small width Suj 
and g(u>) is the density of states defined by Eq. ©. 
Then, we obtain from Eq. ([2^)1 : 

oo 

(e<n in \t)) =\f duuj 2 g{u) lim {N{{Q^Q^f)) . (27) 





Since e„ ^(t) = e„^(i), for iV = oo and all realizations of {-fCre} the time and disorder averaged energy profile is 
given by 



(e n (t)} = 2(ei kin) (0) , (28) 

i.e. the calculation of (e„(t)) is reduced to that of g(ui) and the "quadratic" correlation function ((Q^ QI^ 1 ) 2 ) for 
N — >• oo. {Qn } is the solution of the eigenvalue equation (JSJ with ^ replaced by uj. 

Before we come to the evaluation of the "quadratic" correlation function, let us return to Eq. ([23|) . Making again 
use of the self-averaging of the density of states we obtain for the second term on its r.h.s.: 

oc 

-- f duLJ 2 g(uj) lim (N((Q^Q^) 2 ))cos2iot . 

4 J N— yoo 


Below, it will be shown that lim {N((Q^Qi^) 2 )) is a finite and smooth function of uj. Therefore, the disorder 

N—¥oo 

averaged second term will converge to zero, for t —> oo, due to g(tj) —> go = const., for uj — >• 0. The same property 
should hold for the disorder average of the square bracket term in Eq. (|25|) . With the density of states g(u>,ui') giving 
the joint distribution for two eigcnfrcqucncies the disorder averaged square bracket term becomes a double integral 

over uj and uj' . Although we do not have a rigorous proof, lim (^((Q^Qlf XQ&^Qfe ))) which is part of the 

N— yoo 

integrand should be a finite and smooth function of uj and uj' . Then, taking the limit N — > oo first, the square bracket 
term should converge to zero for t — > oo. If this is true the disorder averaged energy profile converges to an asymptotic 
profile for t — > oo which is consistent with our numerical result. Indeed, the disorder averaged profile in Fig. Q] depends 
on t only very weakly, for large t. In that case the asymptotic profile equals the time averaged one. 

Now we come back to the "quadratic" correlation function. Due to the disorder average it will depend only on 
\n — no|. Since the Anderson modes are exponentially localized one expects that this correlation function decays 
exponentially with \n — uq\. To prove this we first present a crude heuristic approach by assuming 

Q0-)«AUxp(-^-^) (29) 

where the "center of mass" of the Anderson mode v is at n v , which is a random variable, depending on {K n }. J\f v is 
a normalization constant. It should be remarked, that the envelope of an Anderson mode Qn^ decays exponentially, 
but not Q { n ] itself. Therefore Eq. ^ is a crude approximation neglecting sign changes of Qn with n. Substituting 
from Eq. (|2"9"]) into the "quadratic" correlation function and using : 
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1 N 



AT 



we get : 



((Qi^Qi^) 2 ) « ^<[coth|- + |n-no|]exp[-|-|n-no|] (31) 
i.e. the "quadratic" correlation function decays exponentially. 

For an analytical calculation of the "quadratic" correlation function in Eq. (f2"Tj) one can use the approach presented 

in Refs. [22T - |24l | . These authors prove that the computation of the correlation functions (QnQno) ancl (\Qn ||Qno I) 
for \n — no | — > oo is reduced to the solution of an eigenvalue problem for an integral kernel. As a result, these 
correlation functions decay exponentially for large | n — uq | with an inverse localization length given by — In | X max (w) | . 
|A ma x(w)| is the largest absolute value of the eigenvalues of the kernel. It is smaller than one. Applying that approach 
it follows for | n — tiq | — > oo 

((Q i n ) Qn > ) 2 ) = «Mexp(-|n-n |/6M) (32) 

with a correlation length ^(w). We note that the eigenvalue problem in form of an integral equation can only be 
used to calculate correlation functions of the Anderson modes and not directly to compute the energy profile itself. 
But the former is needed (see Eq. (|27p ) for the latter. 

The correlation lengths (localization lengths) of the correlation functions calculated in Refs. (22l42^ | and of the 
"quadratic" correlation function Eq. (f3"2"j) arc different from each other and different from £(w) (Eq. ([TTT0 . for finite 
u>. But for uj — > they exhibit the same divergence, i.e. it is (see Eq. (jlip ): 

f 2 (w) = c 2 uj- 2 , lu^O (33) 

with a positive constant C2, depending on p(K). 

The pre-exponential factor a(u) can be determined as follows. Assuming that Eq. (|32p is valid for all |n — no\, 
summation of the l.h.s. and r.h.s. of that equation and accounting for the normalization XXQ™^) 2 = 1 f° r CJ = 

n 

(remember that Q„ has been chosen as real) yields for N — >• oo: 

a{uj) = 



iVcoth(l/6(w)) 

a |1 , .-,0 . (34) 
In the last line, Eq. Q has been applied. With Eqs. ([52]), Qgj} and (gTJ), it follows from Eq. (ggj: 

OO 

/ — 77T\ - 1 /"^ c ^ 2 ex P(-|" - n o|/6(w)) 

(e ra (t)) = - J du 9 { U ) U CQth(1/6H) • (35) 



The asymptotic \n — no|-dependence is governed by the small-w behavior of the integrand. From Eq. (JTUJ) we get 

g(uj) = dl(cj)/duj ^(K-^/ir , uj -> . (36) 



Assuming that Eqs. (j3"3"]) and (|3l))) are valid for all oj will not influence the asymptotic dependence of (e n (t)) on 
|n — no|. Then we get from Eq. (|35|) for the infinite chain and a displacement excitation: 

CT)> = ^V^" 1 }/* |n-n |- 5 / 2 , |n-no|^oo, (37) 
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i.e. the time and disorder averaged energy profile decays as a power law in \n — no\ with an exponent r\ = 5/2. 
So far we have discussed the energy profile for a displacement excitation. The corresponding calculation for a 
momentum excitation is similar. With the initial condition Eq. (fl~6|) and B = 1, Eq. ([7])leads to 



JV-l 

%W=E Wo )cos ^ " ( 38 ) 

Besides cosw„7j the main difference to u n (t) for a displacement excitation (see Eq. (|2"3")l ) is the absence of the 
prefactor uj v of QnQno- As a consequence one obtains 



/ 7TT\ _ 1 f, 1 ^ cxp(-|7i-n |/^ 2 (^)) , Q , 

(e»(t)) = - J <L,g{ u ) CQth{1/Uu)) (39) 



where ui 2 in Eq. (|35[) is replaced by one. With the same assumptions as above we obtain for the infinite chain and 
a momentum excitation: 



- (K~ 1 )/tt |n-n |- 3 / 2 , |n-n |^oo . (40) 

o 

It is not surprising that we find a power law decay again. The corresponding exponent is rj = 3/2. 

Figures [2] and [3] report the numerical result for the disorder averaged energy profile at different large times of a 
displacement and momentum excitation, respectively. They clearly demonstrate first that, the numerical result of the 
disorder averaged energy profile becomes independent of t for t large enough, and second, that it converges to a power 
law for large \n — uq\ with exponents predicted by the analytical calculation. The three spikes at the i-dependent 
positions n(t) are the phonon fronts propagating with the effective sound velocity, Eq. (|12[) . 




FIG. 2: (Color online) Energy profile at three different times averaged over 10 3 realizations of the disorder with R = 4 for 
iV = 8192 particles and a displacent excitation with A = 2. The dashed line is the predicted power law decay, Eq. (|37|) . The 
ballistic peaks propagate at a velocity c = 1.476 in agreement with the value c = 1.471... computed from Eq. (|12|) . 



There is a finite size effect for N < 00, due to the existence of extended states. For to < uj^ = 0/y/N, (with 6 being 
a suitable constant 0(1)), it is: 

{{Q^Qnof) ^^m 2 (cqn) S m 2 (cqn ) (41) 

where u) = cq. With g{oS) ~ go for the density of extended states it is easy to estimate the contribution of those to 
(e n (t)) in case of a displacement excitation: 



M*)> 



(oxt) 



'5o 



AT-5/2 



(42) 
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FIG. 3: (Color online) Same as Fig. [2] but for a momentum excitation with B 
dashed line is the predicted power law, Eq. (14011 . 



2 and averaged over 2 x 10 realizations. The 



If 1 < \n - no | < N, then it is (e n (i)) = (e„(i)) (loc) + 
but \n — uq\ = 0{N) there is a crossover value \n — no\ c .o. 



3n(«)> 



(ext) r~j 



2n(*)> 



(loc) 



71 - 7l 



" 5 / 2 . For |n - n \ < N 



depending on 9, go etc. such that 



(e„(t)) (loc) = (e„W) (cxt) ^iV- 5 / 2 



(43) 



For N = 8192 this contribution is of order 10 . The corresponding contribution for a momentum excitation is 



(e„(t)) (cxt) - 7V~ 3/2 



(44) 



being of order 10" 6 for N = 8192. 



Figure |4] compares the energy profiles for different strengths of the disorder, i.e. for various values of the parameter 
R. We limited ourselves to the case of a displacement excitation. The profiles display the same decay law. The cases 
with stronger disorder attain the asymptotic profile at smaller distances since in this case the localization lengths are 
shorter. As seen from Figure HJ the data are consistent with the expectation that the asymptotic profile is reached 
for |n — no | 3> £mm- The values of £ m in given in that figure are a rough estimate of the shortest localization length 
obtained by extrapolating the formula pip at uj = uj max = 2c, i.e. £, m in — Ci^max)- 




FIG. 4: (Color online) Disorder averaged energy profiles at t = 2000 for a displacement excitation with A = 2 and increasing 
disorder strengths (top to bottom). Other parameters as in Fig[T] The dashed line indicates the predicted power law Eq. (|37[) 



The prefactor of both power laws, Eqs. (|3"T|) and (|4T)|) depends on the disorder as demonstrated by Figure 4. It seems 
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reasonable that £2 — A£(w) for uj — > with a positive parameter A, independent on u> and the disorder. Eqs. (jTTJ) and 
(f33l) . together with this hypothesis, imply: 



= 8A 



(K- 



(K- 2 ) - (K- 1 )- 



Again m = 1 and a = 1 has been used. For the uniform distribution p(K), Eq. (fl7)) , (if" 
be calculated. From this we obtain: 



(45) 

and (K~ 2 ) can easily 



c 2 {R) = 8 A 



(R- l)RlnR 
(R-l) 2 - R(\nR) 2 



72\(R-1)- 2 [1 + 0(R-1)} 



(46) 



Note, that C2(R) diverges in the no-disorder-limit R — >• 1 + , as it should be since only extended states exist thereby. 
Accordingly, should become infinite for all u>. 

Introducing c 2 (R) from the first line of Eq. (ggj and (if _1 )(J?) = (In i?)/(i? - 1) into the prefactor C(i?) of the 
power laws, Eqs. (|57|) and (|4U)) . leads to the independence shown in Fig. [5] in case of a displacement and a momentum 
excitation, respectively. The unknown parameter A has been adjusted in order to fit the numerical result for the 
prefactors. The latter are obtained from the numerical data in Figures [H and [3] extrapolating (e„(t)) \n — no\ n at large 
\n - n |. 



Displacement exc. 

(\ lt = 3 - 97 > 
Momentum exc. 
= 144) 




R 
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FIG. 5: (Color online) Dependence of the prefactors C(R) on the disorder strength R for a displacement (A — 2) and momentum 
(B = 2) excitation: numerical (full circles and triangles respectively) and analytical result Eqs. (|37p . (|4Up with formula (|46p 
(solid lines). Dash-dotted lines are the expected asymptotic behaviors for R — > 1 + , (R— l) -3 and (ii— 1) _1 , respectively. The 
fitting parameter Aflt is given in the legend. 

The numerical and analytical result for the prefactor in case of a displacement excitation agree satisfactorily, even 
for the smallest value of R = 1.3. Investigating the profile for even smaller values is hampered for our finite chain 
by the increase of the localization length with decreasing (R — 1). The same agreement is also valid in case of the 
momentum excitation, except for the two smallest i?-values at 1.3 and 1.5. Eq. ([7]) demonstrates that the weight of the 
low- lying Anderson modes for a momentum excitation is by a factor \/u) v higher than for a displacement excitation. 
Since the localization length increases with decreasing uj v , this could be the reason for the "asymmetric" behavior of 
C(R) for both kind of excitations. Indeed, we have observed that (e n (t))\n — n a \ z / 2 does not reach a stationary value, 
for e.g. R = 1.3. The strong deviation of the fit parameter A for the displacement and momentum excitation may 
originate also from this fact. 



C. Moments of the local energy 

A customary way to describe wavepacket diffusion is to look at time evolution of moments of the energy distribution 
that are defined as 

m v {t) = Snl"-"oren(*) (4?) 
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(the denominator is clearly only a scale factor) . Of particular interest for a statistical characterization are the disorder 
averaged moments (m u {t)). Their numerical result is shown in Figure [51 If one uses the asymptotics Eqs. ([57)1 , P0|) 
and introduces a cutoff of the sum in the numerator of Eq. (|47p at the ballistic distance \n — hq\ = ct one obtains 

(m„(t)) oc ftu) = { V Q + 1 - f] \ V v >lZ\ . (48) 

For v = rj — 1 there is a logarithmic divergence of (m v (t)) with time. As demonstrated in Fig. the numerical 
values of (3{v) are in excellent agreement with Eq. (|48|) . This also implies that the contribution of the traveling peaks 
is not relevant as implicitly assumed in the derivation of Eq. (|48j) . 

This result is consistent with the values that could be inferred by Datta and Kundu [l(|. Indeed, they predict 
/3(2) = 1/2 and /3(2) = 3/2, respectively, for a displacement and momentum excitation. Notice that, if one looks 
only at m,2 (t) one would incorrectly conclude that the two cases would correspond to sub- and superdiffusive behavior 
respectively. A full analysis of the spectrum of moments and of the wavefront shape is necessary to assess the real 
nature of dynamics. 




time 

FIG. 6: (Color online) Evolution of moments {m v {t)) for v — 0.5, 1, 1.5, 2, 2.5, 3 (bottom to top respectively) for the harmonic 
chain TV = 8192, R — 6. Displacement (a) and momentum (b) excitation (A = 2 and B = 2 respectively). Averages are over 
6 x 10 3 disorder realizations. 



IV. ENERGY PROFILE: ANHARMONIC CASE 

In this section we will investigate numerically the n-dependence of the energy profile averaged over the disorder in 
the presence of anharmonicity. Particularly, we will check whether its tails can be described by those of the harmonic 
chain. As a model we have chosen the Fermi-Pasta-Ulam (FPU) chain with cubic nonlinear force 

mu n = K n (u n+1 - u n ) - -fT n _i(it„ -u n -i) 

+G(u n+ i - u n f - G(u n -! - u n ) 3 . (49) 

It reduces to the harmonic chain for G = 0. For simplicity, we considered the case of uniform nonlinear coupling G 
{G = 1 in the following). 

The analysis of the previous section shows that the behavior of the harmonic chain follows all the expected features. 
Which influence of the anharmonicity do we expect? If the initially localized energy would spread completely it would 
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FIG. 7: (Color online) Comparison between the exponents measured numerically for momentum and displacement excita- 
tions (full circles and squares respectively) and the analytical result Eq. (|48|l , same parameters as in Fig. [6] The exponents 
are evaluated by a power-law fit. Error bars are estimated from the fluctuations of the (discrete) logarithmic derivative 
A log{m„(f)) /A log t and are reported only when larger than symbols' size. 

be e n (t), — !► for t — > oo, for all n. For incomplete spreading, however, (e„(i)) for t large enough should decay by 
the power laws Eqs. ([3"T]) or Eqs. (|40p . again, and the amplitudes of oscillations at sites far away from site no of the 
initial excitation should become so small that the harmonic approximation applied to those tails should become valid. 

A detailed analysis of the effects of nonlinearity goes beyond the scope of the present work. We thus limited 
ourselves to the case of FPU with initial displacement excitation with A ~ 2. Wc checked that the energy is about 
a factor of 2 larger with respect to the G = case meaning that the nonlinear part of the potential is sizeable. Wc 
considered the usual definition of e^ *-* where Ki[ui+i — itj], i = n — 1, n in Eq. (fT4")l is replaced by V^(ui+i — u{) with 
Vi{x) = KiX 2 /2 + Gx A /4. As for the harmonic case, we performed the average over disorder at three different times. 

The average energy profiles for three different disorder strengths are reported in Fig. [5J The profiles still show a 
pretty slow decay, reminiscent of the harmonic case. From Fig. [SJwe first observe that the convergence of (e n (t)) to a 
limiting profile at t = oo becomes slower for larger disorder strength R, i.e. for shorter localization lengths. Second, 
whereas the profile for R — 4 and the largest time t = 6000 can be satisfactorily fitted by the power law Eq. (|5T| , this 
is less obvious for R = 2 and R = 8. For R = 2 and \n — Uo\ > 1000 the profile is practically time independent for 
t > 2000. But in contrast to the harmonic case (see Fig. [2]) it does not reveal the asymptotic power law \n — n \~ 5 / 2 , 
although the data suggest that this may happen for \n — n$\ > 3000. For R = 8 the profile follows that power law 
for 100 < \n — no| < 1000, i.e. for about a decade, but deviates for \n — tiq\ > 1000. However, comparing this profile 
for the three different values for t hints that the range of the power law decay may increase with increasing t. In 
addition, the profiles display some form of weak " broadening" of the tails indicating that some energy is indeed slowly 
propagating. 

As a consequence, the disorder averaged moments (m u {t)) do not display a convincing scaling with time. Even for 
statistically accurate data as the one in Fig. [51 the effective exponents (as measured for example by the logarithmic 
derivatives of (m„(i))) display sizeable oscillations which are well outside the range of the statistical fluctuations (see 
Fig. [!)]). Similar results are obtained for momenta of different order (not reported). 

We may thus argue that, at least in the considered parameter range, the nonlinear case has a core which remain 
almost localized (in a similar way as the harmonic case) but in addition there must be a small propagating com- 
ponent. The fraction of such propagating component increases upon increasing the energy and/or nonlinearity. As 
a consequence, with the data at hand it is impossible to draw definite conclusions on the nature of the spreading 
process. 

V. SUMMARY AND CONCLUSIONS 

The relaxation of an initially localized excitation in a translationally invariant chain of particles has been studied 
for harmonic and anharmonic nearest neighbor couplings. The main focus has been on the energy profile (e„(i)), the 
moments (m^(t)), both averaged over the disorder, and the relation between the asymptotic t-dcpcndcncc of (m u (t)} 
with the asymptotic n-dependence of (e n (t)}. As far as we know, this has neither been explored for the anharmonic 
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FIG. 8: (Color online) FPU model: Disordered-average energy profile at three different times averaged over 10 realizations 
of the disorder and for different disorder strengths R. Chain of N = 8192 particles with displacement excitation A = 2. For 
comparison, the predicted power law decay for the harmonic chain, Eq. (|37[) . is also drawn (dashed lines). 




FIG. 9: (Color online) FPU model: Evolution of (m2(i)) for R = 2, 4, 8 (top to bottom respectively) same parameters as in 
Fig. [8] The inset reports the (discrete) logarithmic derivative A log(m2(t))/A log t of the data versus logi. For comparison, 
the value for the harmonic chain /3(2) = 1/2 is also drawn (dashed horizontal line). 



nor for the harmonic case due to the lack of analytical knowledge of (e n (t)) for |n — uq\ — > oo. 

For the harmonic model we succeeded to determine analytically the disorder and time averaged energy profile (e n (t)) 
for a displacement and a momentum excitation at site no and initial time t — 0. Whereas e n (t) is a quasiperiodic 
function which does not converge for t — > oo we have argued that (e n (t)) converges for t — > oo. In that case (e n {t)) 
gives the limiting profile averaged over the disorder. The analytical calculation yields a power law decay 



(e n {t))^C{R)\n-n \-* 



(50) 
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for 1 <C \n — no | <C N, in case the system is finite. The exponent rj and the prefactor C(R) depend on the type 
of excitation. For a displacement and momentum excitation we have found rj = 5/2 and rj = 3/2, respectively, in 
good agreement with the numerical values. This agreement also holds for the analytical and numerical results for the 
i?-dependence of C(R), except for the two smallest values of R in case of a momentum excitation. Accordingly our 
assumption ^(w) ~ for uj — > is supported. From this proportionality it also follows that C(R) diverges at 
R = 1, the no-disorder limit. The power law decay, Eq. ([50j). originates from the gapless excitation spectrum of the 
Anderson modes. It is the consequence of the translational invariance of model ^ . Destruction of this invariance 
by adding, e.g. an on-site potential like in the KG model generates an energy gap. The corresponding localization 
length at the lowest eigcnfrcqucncy will not diverge anymore, and therefore the energy profile will decay exponentially 
for | n — KqI — > oo. However, we stress that any lattice model without an external potential has to be invariant under 
arbitrary translations. This implies a gapless spectrum which is the origin of the power law decay of the profile. 

The power law decay of (e n (t) has remarkable consequences on the asymptotic ^-dependence of the moments. If 
we use Eq. (|50[) to calculate the time and disorder averaged i/-th moment we get for displacement and momentum 
excitation: 



{m v (t)) = oo (51) 

for all v > 2; note that, for instance, (mi(t)) is finite for a displacement, but not for a momentum excitation. The 
result, Eq. (|5ip. implies that the disorder averaged moment (m v (t)} must diverge with time, although the initial local 
energy excitation does not spread completely. This power law divergence of (m„(i)) with time is clearly supported by 
the numerical result for v > 2 and v > 1 for the displacement and momentum excitation, respectively. As a matter 
of fact, consideration of m^it) alone is not sufficient to conclude that the energy diffuses. This is one of the main 
messages of the paper. 

The analytically exact result for (e n (t)) in case of harmonic interactions also allows to check how far the tails of 
an anharmonic chain, where the average displacements become arbitrary small, can be described by the tails of the 
harmonic system Although no definite conclusion can be drawn, we have found evidence for a crossover of the energy 
profile of the anharmonic to that of the harmonic chain. However, for the weakest and strongest strength of disorder 
this crossover seems to occur for \n — n \ > 3000 and for t > 6000, respectively. This may be explained as follows. 
The localization length ^(w) is large for weak disorder. Since \n — no|/^2(w) enters into the calculation of the disorder 
averaged profile (see Eq. (|35[0 the asymptotic power laws, Eqs. (|3"T|) and (|4T)|) . occur at larger values of \n — Uq\. 
For large disorder, ^(w) is small. But the time scale for tunneling processes responsible for the energy propagation 
increases significantly due to an increase of the potential barriers. Therefore the convergence to a limiting profile is 
much slower which is exactly what we observed (see Fig. [7]). The increase of the localization length for weak disorder 
and the increase of the relevant time scale of the relaxation for strong disorder probably are also the reasons for the 
absence of a convincing scaling of the moments (m v (t)). In order to test this, one has to increase both, the number 
of particles and the simulation time significantly. Requiring a similar good statistic of the data this has not been 
possible so far within the available CPU time. If it is true that the asymptotic energy profile agrees with that of the 
harmonic chain this would imply that the moments (m v (t)) for the anharmonic model for v > 2 diverge with time, 
as well, although the energy does not spread completely. 

From our results, it is nonetheless clear that the interplay of localized and almost-extended modes leads to a 
nontrivial decay of wavepackets amplitudes and this must be taken into account when dealing with the nonlinear case. 
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